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Abstract 

Electromagnetic material characterization is the process of determining the con¬ 
stitutive parameters of matter. In simple media, these parameters are permittivity 
and permeability. Characterization of these values is often accomplished through the 
use of waveguides, transmission lines, coaxial cables, and resonant cavities. Free space 
measurements systems are employed since they are non destructive (i.e., no sample 
machining is required) and broadband. Traditionally, time domain gating is utilized 
to mitigate systematic errors. However, an artifact of this calibration technique is 
band edge corruption due to data windowing. The goal of this research is to de¬ 
velop and apply a Three Short Calibration Technique to the General Electric Low 
Frequency (0.5-2GHz) Focus Beam System in order to preserve band edge data. It 
is shown that when working with low frequencies, such as the focus of this research, 
coupling effects due to multiple bounces within the collimating and focusing lenses as 
well as lens-sample interaction are not easily calibrated out. Although the coupling 
effects introduce extraneous clutter into the extracted constitutive parameter values, 
the results obtained agree with the known reference values of tested dielectric sam¬ 
ples. It is concluded that with a more in depth study in sample/lens interaction at low 
frequencies, this calibration technique could prove useful in the accurate extraction 
of the desired low frequency constitutive parameter values. 
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APPLICATION OF THE THREE SHORT CALIBRATION TECHNIQUE IN A 


LOW FREQUENCY FOCUS BEAM SYSTEM 

I. Introduction 

The measurement of a sample’s constitutive parameters is of great interest to 
engineers as it allows for the accurate prediction of how electromagnetic waves will 
interact with a particular material. Evaluation of the complex permittivity and per¬ 
meability is required not only for engineering but also industrial applications. Nu¬ 
merous material characterization techniques exist, each with particular advantages 
and disadvantages. There is no single optimum technique for all material charac¬ 
terizations. In trying to characterize low observable materials, traditional methods 
such as waveguides, transmission lines, co-axial cables, and resonant cavities tend 
to introduce error due to inhomogeneities, air gaps, and finite conductivity that are 
difficult to correct. Free-space techniques, such as the focus beam system, offer an 
advantage for characterization when dealing with these types of materials since they 
are broadband and nondestructive. 

1.1 Problem Statement 

An automated network analyzer system is used to make frequency swept S- 
parameter measurements on the focused-beam sample fixture from which complex 
permittivity and permeability may be calculated. However, before any measurements 
can be taken the NWA must first be calibrated. 

Automated network analyzers inherently have many errors associated with the 
measured scattering parameters. The sources for these errors according to [1] are 

• Imperfect connector matching 

• Imperfect calibration standards 

• Nonlinearity of mixers 
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• Gain and phase drifts in IF amplifiers 

• Noise introduced by the analog to digital converter 

• Imperfect tracking in dual-channel systems 

Calibrating the NWA attempts to remove the systematic uncertainties by mea¬ 
suring a set of known standards. After calibration, the NWA operates with error 
correction and automatically updates the measurements with the calibration infor¬ 
mation. The complex material parameters can then be computed from the measured 
scattering parameters using a computational algorithm such as the Nicolson-Ross- 
Weir (NRW) or Newton-Raphson Root Search techniques [12,20]. 

The Low Frequency Focus Beam System at General Electric currently uses a 
time domain gating technique in order to obtain calibrated S-Parameter measurements 
which subsequently leads to the extraction of material constitutive parameters. The 
problem with this method is that critical data is attenuated in the high and low 
frequency region (i.e. band edge data) due to data windowing. General Electric is 
interested in preserving the data that is lost through this process, especially at the 
low end of the frequency band. The three short calibration technique aims to be an 
effective way of fully calibrating this system, so that the time domain gating method 
may be avoided [18]. 

1.2 Limitations 

The three short calibration technique is band limited due to degeneracy of shorts 
spaced greater than or equal to a half wavelength. Special care must be taken in 
the derivation of the maximum spacing constraint between the shorts. This will be 
accomplished by establishing a set bandwidth for the experiment (0.5-2GHz). 

General Electric uses a mechanical mechanism to reposition the sample holder. 
Although this is an extremely good idea, accurate precision may prove to be an issue as 
a high degree of accuracy when translating from position to position is needed in order 
to not introduce phase errors in the extracted data. It will also be important to see 
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if this change in phase can be resolvable. There are also other errors associated with 
measurements in the low frequency regime, including lower signal levels embedded 
in low frequency clutter, coupling effects, and antenna/lens internal reflections. The 
three short calibration technique investigated in this research seeks to mitigate these 
errors. 

1.3 Scope 

The three short calibration technique will be applied to the Low Frequency Focus 
Beam System. As previously discussed, this calibration technique is band limited. 
Communication with the engineers at General Electric established an interest in the 
bandwidth between 500Mhz and 2GHz. Two dielectrics will be tested for validation, 
a fiberglass panel and a sheet of plcxiglas. The current technique of time domain 
gating will be used on the fiberglass panel to establish a comparison of data with the 
full frequency domain calibration results. An analysis of post calibration error will be 
discussed. 

1.4 Thesis Organization 

This thesis is divided into five chapters. Chapter two provides a brief description 
of permittivity and permeability. Current calibration techniques used on similar Free 
Space measurement systems will be discussed and pros and cons will be reviewed. A 
section on the time domain gating technique will also be reviewed to increase knowl¬ 
edge in the need for an alternative calibration technique for this system. Chapter 
three provides an overview of the NWA 12-term error model, since the 3-short cali¬ 
bration technique is used to calculate these error terms. The maximum short-spacing 
constraint is established and a verification of phase angle resolution is performed. 
Chapter four discusses the results of the experimental measurements as well as specu¬ 
lation on any error still present after application of the calibration technique. Chapter 
five will give a conclusion to all results found in this thesis. It will also suggest changes 
to the experimental setup that can be made to possibly enhance accuracy. 
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II. Literature Review and Background Material 

I n order to understand the importance of measuring the constitutive parameters 
of materials, or even the reason behind taking an interest in developing multiple 
techniques in doing so, one must strive to develop a fundamental understanding in 
certain key areas involved in the held of electromagnetics. This chapter will focus on 
background concepts that will bring out theory and solidify the justification of why 
the ” three short calibration technique” can be applicable to the focus beam system. 

2.1 Properties of Matter 

Materials are typically defined and described by their constitutive parameters. 
The main purpose of most material measurements is to ultimately obtain data that 
will allow for calculation of the material’s permittivity and/or permeability values. 
These values are known as a material’s constitutive parameters. There is no method 
to date that allows for the direct measurement of a material’s constitutive parameters, 
these values are calculated through various novel numerical and iterative techniques, 
each having their own advantages and disadvantages. Once a material’s constitutive 
parameters have been extracted, fairly accurate prediction of how the material will 
react when an electromagnetic wave becomes incident upon it can be determined. 
This information is of importance in the engineering community because it allows for 
engineers to place layers upon perfect electrical conductors (PECs), for example, in 
an effort to control electromagnetic scattering. 

A further discussion of permittivity and permeability is needed in order to un¬ 
derstand the big picture of going through the trouble of extracting these values. In 
order to develop these definitions a discussion on the physical structure of matter will 
be presented. A review of these topics will enhance the understanding of the vari¬ 
ous methods for extracting the constitutive values of materials, and ultimately the 
understanding of the importance of advancing the calibration procedure for the Low 
Frequency Focus Beam System. This thesis involves the application of a calibration 
technique on a system that will be utilized to measure the scattering matrix of a 
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E=0 



Figure 2.1: Polarization of a nonpolar atom or molecule. 

device under test. This scattering matrix will in turn be mathematically manipulated 
in order to extract permittivity and permeability values of the sample of interest. 

2.2 Dielectrics 

Following a development in [21], a concise derivation of the dielectric constant 
is presented. In order to break down the structure of a dielectric, its atom will be 
depicted by a negative charge — Q (electron cloud) and a positive charge +Q (nucleus). 
The nuclei will represent point charges and the electron cloud consisting of multiple 
negative charges. There is assumed to be an equal amount of positive and negative 
charges, making the atom electrically neutral. When an electric field E is incident 
upon this atom, the negative charge is displaced from its equilibrium position as 
shown in Figure 2.1. The result of these actions is the creation of a dipole, and the 
dielectric becomes polarized with dipole moment 

p = Qr (2.1) 

where r is the distance vector from — Q to +Q. This is depicted in figure 2.1. If there 
are N dipoles within a volume Sv of the dielectric, the total dipole moment due to the 
electric held becomes 

N 

Qi'E + Q 2 r 2 + ••• + QnTn = QkJk- (2.2) 

k =1 
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The density of these dipole moments therefore can be expressed as 


P 


lim 

Av—^0 


N 

^ Qk^k 


k =1 


Av 


(2.3) 


The major effect of E on a dielectric is the creation of dipole moments in the direction 
of E. The behavior previously described is that of a nonpolar dielectric. Dielectrics 
with built in permanent dipoles that are randomly oriented are called polar dielectrics. 
When these types of dielectrics are in the presence of E, their dipoles align themselves 
with E. leading to a net polarization density of p. 

When applying E to a dielectric material, the electric flux density D is greater 
than it would be in free space. This is defined by the relation 


D = e 0 E + P. 


(2.4) 


If a. material is linear and isotropic, the polarization will be proportional to the applied 
electric field, thus 

p = XetoE (2.5) 

where \ e is the electric susceptibility of the material. By substituting equation 2.5 
into equation 2.4, one can now define D to be 

D = eE (2.6) 


where 

e = £oe r (2-7) 

and 

e r = 1 + Xe = —• ( 2 - 8 ) 

eo 
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Figure 2.2: Electron orbiting around the nucleus 

The parameter e is said to be the permittivity of the dielectric, eo is the permittivity 
of free space having value 

p 

e 0 « 8.854 x 10” 12 —, (2.9) 

m 

and e r is the dielectric constant or relative permittivity. 


2.3 Magnetics 


A similar development presented in [21] will be followed for the description of 
relative permeability. In parallel with the development of polarization, electrons also 
rotate about their own nucleus. The act of an electron orbiting around the nucleus 
creates a magnetic dipole m as illustrated in figure 2.2. 


Without the presence of an external B field, the sum of magnetic moments m 
for many materials is zero due to random orientation. This phenomena is depicted in 
Figure 2.3. In the presence of an external B field, the magnetic moments leading to 
the magnetization density 


M 


lim 

A-u—>0 


N 

E 


k= 1 


Av 


( 2 . 10 ) 
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Figure 2.3: Magnetic dipole moment before and after B is applied. 

Similar to D, an auxilliary fied H is defined to accommodate the effects of the mag¬ 
netization density, namely, 

H=— ( 2 . 11 ) 

Vo 

or 

B — Vo{H + Vo) ( 2 - 12 ) 

, where H is the magnetic field intensity. 

Equation 2.12 holds for all materials whether they are linear or not. For linear 
materials, M depends linearly on H that is 

M = XmH (2.13) 

where Xm is the magnetic susceptibility of the medium. Substituting equation 2.13 
into equation 2.12 leads to 


B - /io(l + Xm)H - fiH 


(2.14) 










or 

B = noHrH (2-15) 

where 

/i r = 1 + Xm = —• (2.16) 

ho 

The constant fi is called the permeability of the material having unites of henrys/me- 
ter. The dimensionless quantity /i r is the ratio of the permeability of a given material 
to that of free space and is known as the relative permeability of the material. 

2.f The Low Frequency Focus Beam System 

In order to further reveal the problem of band edge corruption, a brief descrip¬ 
tion of the system to be evaluated is needed. The system under evaluation is known 
as the Low Frequency Focus Beam (LFFB) System. In the context of material charac¬ 
terization, this particular system was designed to emulate plane-wave measurements 
of the transmission and reflection coefficients from planar samples where the solutions 
for intrinsic properties can be easily derived [5,6,19]. 

The LFFB system consists of horn antennas, collimating/focusing lenses and a 
sample holder, as depicted in Figure 2.4. The planar sheet to be measured is placed 
in the focused beam, typically at the beam waist of focus. This minimizes scattering 
do to the illumination of material edges, which is a well known inherent source of 
error for this particular type of system. The beam waist for this system varies as a 
function of frequency. Note, the reflection and transmission coefficients, which are 
vital for extracting e and //, are measured using a network analyzer. 

2.5 Types of Calibration Techniques 

The data that is recorded by the vector network analyzer is collected at the 
internally-housed detector elements. In order to acquire the sample/material scat¬ 
tering parameters, this data must be calibrated in order to account for the inherent 
systematic errors that are associated with the instrumentation and test fixture. Many 
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Figure 2.4: The General Electric Low Frequency Focus Beam System 

calibration techniques exist, each with their own obvious advantages and disadvan¬ 
tages. This section will focus on some of the more popular calibration techniques that 
have been used. The bulk of calibration processes can be separated into two groups. 
These groups are known as Standard Calibration and Self Calibration, as discussed 
next. 


2.5.1 Standard Calibration. The Standard Calibration is commonly known 
as the Short-Open-Load-Through (SOLT) technique. The downside to this calibra¬ 
tion technique is that a large number of standards have to be measured, which can be 
time consuming. In addition, these standards may be difficult to achieve in practice. 
For example, when applying this calibration technique, as the frequency increases, 
the short and open standards are difficult to define because of the increasing effect 
of parasitic inductances and capacitances [3]. The SOLT calibration method is also 
rarely used for Free Space measurements because the Open standard is difficult to 
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fabricate, although attempts have been made using artificially-created periodic struc¬ 
ture technology. Thus, for free-space or focus-beam systems, alternative calibration 
techniques are required, as described next. 

2.5.2 Self Calibration. Most other calibration methods fall into the category 
of Self Calibration techniques. These types of calibration techniques are different from 
Standard Calibration techniques because they use data obtained from transmission 
lines rather than data gained from discrete standards [10]. Self Calibration techniques 
can be broken down into two categories, Through-Reflect-Line (TRL) calibration and 
Through-Reflect-Match (TRM) calibration. 

The TRL calibration technique utilizes a transmission Line of known length 
and impedance as one standard. This Line standard has to be significantly longer in 
electrical length than the Through line. For radio frequency applications, the length 
of transmission lines needed to encompass low frequencies become impractical [22], 
TRL calibration also involves the use of a high-reflection standard which must be 
electrically the same for both test ports. Although this type of calibration technique 
when applied to the free space measurement system can result in error clue to move¬ 
ment of actual components in the system, this method has been applied successfully 
in various experiments [6,15-17]. Gagnon presents a study on the Misalignments in 
Quasi-Optical Measurements [4], In reviewing this, it is concluded that a calibration 
technique that does not involve translation of the feed horns or movement of the sam¬ 
ple lens is strongly advised. Slight misalignments cause error to manifest itself within 
the results that are hard to correct. Previous experiments suggest movement of the 
feed horn in order to realize the “Line” standard. A calibration technique that does 
not involve movement of the feed horns is sought after for the purposes of minimizing 
systematic errors. If the line standard is to be used, the difference between the two 
lengths of line must be approximately a quarter of wavelength at midband frequency, 
therefore the TRL is bandwidth limited [10]. In many cases of Free Space Measure¬ 
ment Systems, the horns are mounted to ensure their precise positioning in relation 
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to the collimating lenses. Thus, movement of the horns in order to record data for a 
’’line” standard is not recommended. 

The need for a better approach in calibrating this system led to alternative 
calibration methods. The Through-Reflect-Match (TRM) Calibration with the addi¬ 
tion of some correction terms described in [5] is another validated calibration process. 
This process is fairly similar to the previously discussed TRL calibration method. In¬ 
stead of using two transmission line standards, the TRM calibration utilizes only one 
transmission line based standard. A highly reflective standard is still needed. The 
major notable difference is brought about by the use of a Match standard instead 
of an additional transmission line measurement. The use of a Match standard is fa¬ 
vorable by some scientists and engineers as it allows for fixed location of the horns, 
therefore less error will be introduced into the results by improper positioning of the 
horns [4]. The match standard also is not bandwidth limited. At the frequencies used 
for the purposes of this thesis, a true match standard is hard to achieve as broadband 
absorbing standards are difficult to fabricate 

2.6 Time Domain Gating Technique 

The time-domain gating calibration scheme utilizes a Reflect setup to initialize 
the reference calibration plane, where the VNA records the S™ 1 response in the fre¬ 
quency domain. This uncalibrated measurement contains systematic error caused by 
the constructive and destructive interference from multiple scatterers present in the 
system (i.e. horns, lenses, etc). Applying a Fourier Transform of this response to take 
it into the time-domain reveals the relative locations of the multiple reflections. This 
exemplifies a physically-intuitive response where the dominant reflection is caused by 
the inserted short while the other reflections are caused by inherent imperfections 
in the system. Multiple reflections from the collimating lenses are a source of these 
reflections, as well as the lenses [7]. These system imperfections are easily identified 
based on the time separation between pulses (see Figure 4.1). 
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The unwanted scatterers can be easily removed by applying time-domain gating 
to the response. This effectively removes the undesired responses and leaves only 
the dominant scatter produced by the reflective standard. Transforming this gated 
response back into the frequency domain produces the desired S response of any 
additional clutter. This time-domain gated short measurement can then be utilized in 
conjunction with the following Thru setup to calibrate the Free Space Measurement 
System. 

The Thru utilizes the Sf[ l response of an empty measurement in the frequency 
domain. Again, the data is transformed into the time-domain, gated and transformed 
back into the frequency domain, producing the desired Sff response that can then 
be used to calibrate the system [7]. 

The two calibration measurements, the Reflect and Thru , can now be related 
to the testing sample measurements, via the relations 

srr = E RT sir , (2.17) 

S™ = E TT S™ m . (2.18) 

E R t and ETT represent the forward reflection and transmission tracking error terms 
respectively, discussed in Chapter 3. The transition region scattering parameters can 
be computed from the calibration measurements. Since the gated short and empty 
measurements are related to the transition region scattering parameters via 

Sfl 1 = E RT Sf; = E rt (-1), (2.19) 

S™ 1 = E TT e~^ l % (2.20) 

where 

k 0 = cuyTo/io, (2-21) 

7o = jk o- (2.22) 
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Where l s is the sample width. Therefore, the de-embeded sample scattering parame¬ 
ters can be obtained through substitution of equations(2.19) and (2.20) into equations 
(2.17) and (2.18), leading to 


Qsam _ _ 


Seal ’ 


Ssam _ 


s? 


SSPerto 1 * 


(2.23) 

(2.24) 


The permittivity and permeability can be subsequently determined using the standard 
NRW [20] technique or Newton-Raphson [12] root search algorithm. 


VNAs that have the ability to do real time Fourier Transforms are expensive 
and some labs don’t have the capability. The main reason for researching an alter¬ 
nate approach to the time-domain gating technique is the fact that you essentially 
have to deemphasize through weighting important band-edge data due to windowing. 
Furthermore, the windowing size itself directly affects the extracted data. Depending 
on window size, more data in the lower and higher frequencies may become corrupted 
due to a weighting factor that is applied on the band edges in order to reduce side 
lobes, in an effort to increase SNR. The aim of this thesis focuses on extracting low 
frequency data, making the time-domain gating method unacceptable. Before pro¬ 
ceeding,to the proposed three short calibration technique in chapter 3, the necessary 
background in material extraction is discussed. 


2 .7 Windowing 

In order to understand the need for a better way to calibrate the Low Frequency 
Focus Beam System, a brief discussion on windowing in the frequency domain will be 
presented. As depicted in Figure 2.5, the use of a rectangular window captures all of 
the frequency data of interest. However, the downside to this window is that it only 
allows for a -13dB drop in the first side lobe level. This small decrease in side lobe 
levels can cause for systematic or environmental noises to be an issue when trying to 
extract important data. Note, the main side lobe width is extremely narrow, which 
can be advantageous in regards to resolution. Also depicted, for comparison is the 
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Figure 2.5: The top picture represents the use of a rectangular window and its 

frequency response. The bottom picture represents the use of a hamming wind and 
its frequency response. This figure depicts the information that is lost due to a trade 
off when trying to decrease sidelobe levels. 

use of a hamming window. It is evident that there is data in the high and low end 
of the bandwidths that is heavily weighted. Although the first side lobe is decreased 
dramatically, it comes at the expense of band edge data corruption. The three short 
calibration technique aims to eliminate the use of windowing all together in an effort 
to preserve the apparent data that is lost as depicted in Figure 2.5 as a result of 
windowing techniques. 


15 
























































































































2.8 Nicolson-Ross- Weir Technique 

The Nicolson-Ross-Weir (NRW) technique allows for the extraction of the com¬ 
plex permittivity and permeability of a single layered material [20]. The algorithm 
itself requires two measured parameters: the transmission and reflection coefficients 
from the material under test. Thus, the forward (Sn, S 2 \) and reverse (S 22 , S 2 1 ) 
scattering parameters can be utilized to compute the material parameters. 

The formulation of the NRW technique compares the theoretical and experi¬ 
mentally measured S-parameters. 

Sn y {u, e, //) - S^ p (uj, e, +) = 0, (2.25) 

S t 2 l l y (u,e,n)-S t 2 2 y (u,e,n) = 0 (2.26) 

Where the theoretical S-Parameters can be formulated according to [6] as 

Sl';’ = -Jpl = Su". (2.27) 

SS” = C 0 -Jpl = SS”. (2.28) 

A closed-form solution can be found by solving these nonlinear equations for R and 
P in terms of S^ p and S^ p . Solving these equations leads to the quadratic equation 


R 2 - 2 QR + 1 = 0 


(2.29) 


Q 


{st x l p y - ( s e 2 xp ) 2 +1 
2 s*? 


(2.30) 


which has the solution 

R = Q ± VQ 2 - 1 (2.31) 


Since the magnitude of the interfacial reflection coefficient can never be greater than 
unity for a passive material, the proper root choice is simply the one which satisfies 
the requirement \R\ < 1. 
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With the experimentally measured S-parameters related to R and P, the complex 
material parameters can now be determined. When utilizing a rectangular X-Band 
waveguide operating in the TE10 mode, assuming no gaps exist between the sample 
and perfectly conducting waveguide walls, the relative permittivity and permeability 
can be formulated as 


z-1 Z 7 o 

R — j 2 — . , — hr 

z + 1 Z 0 jup o/ 7 o 7 

In a free space measurement system, k c = 0 , thus 


(2.32) 


7 2 = k 2 — k 2 , k 2 = co 2 ep, k 2 c = 0 (2.33) 

1 71 "P / 

P = e" 7 ', 7 =-— = y 'k 2 c - k 2 e r p r (2.34) 

Therefore, e r and /i r can easily be extracted using the following two equations. 


7 InP, 1 + R , 

ii r = z— = - - -- 

7o 7o^ 1 - -R 


A ; 2 - 7 2 _ _k 2 c - {InP/l) 2 InP 1 + R, 
fi r kl A; ( j 7 0 / 1 — 


(2.35) 

(2.36) 


Newton Raphson Method 

This section focuses on a numerical technique that will be utilized to extract 
the complex permittivity (or permeability) of the device under test. Many numerical 
methods exist to determine the roots of a nonlinear equation. The Newton-Raphson 
method is a derivative-based root finding technique that only requires an initial guess 
of the root [ 12 ], Convergence is not guaranteed, but if the method does converge, 
it does so much faster than in similar numerical bracketing methods, such as the 
bisection method [13]. 

The Newton-Raphson method is based on the principle that if the initial guess 
of the root of f{x ) = 0 is at ay, then by drawing a tangent line to the curve at 
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f(xi), the point Xi + 1 where the tangent crosses the x-axis is an improved estimate of 
the root. The method can be derived by using the Taylor expansion of an arbitrary 
nonlinear function f(x) 


f(x + Ax) = f(x) + f'(x) Ax = 0 


(2.37) 


Solving for Arc leads to 

Ax 


which can then be expanded as 


fix) 

fix)' 


Arc = rcj + i — Xi = Ax 


fjXj) 

rixiY 


This can then be rearranged to produce the following 


%i -\-1 %i 


fjXj) 

fixtY 


(2.38) 


(2.39) 


(2.40) 


The preceding equation is known as the Newton-Raphson formula for finding the roots 
(i.e., zeros) of nonlinear functions. 


The Newton-Raphson Method can be modified to find the dielectric constant of 
a non-magnetic material. This is accomplished by setting the nonlinear function to 
either 


/M 



(2.41) 


where tol is a user defined tolerance (typically 1x10 7 ). If the derivative cannot be 
found analytically, the numerical derivative of /(re) can be approximately computed 


as 


fix) 


lim 

A#—>-0 


/(X + Arc) - /(re) 

Ax 


/(x + Ax) - /( x) 

Ax 


(2.42) 


The theoretical scattering parameters can be found through the use of wave 
transmission matrices, or A-parameters which describe the relationship at particular 
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Figure 2.6: Reflected and transmitted waves at a planar interface 


reference planes between incident and reflected wave amplitudes. The general A- 
parameter formulation considers incident waves C\ and b' 2 on an interface from the left 
and right, respectively, as depicted in Figure 2.6. The respective interfacial reflection 
and transmission coefficient experienced by waves Ci and b 2 are represented as i? 1; T 12 
and R- 2 , T 2l such that 

c 2 = T 2 \Ci + R 2 b 2 (2.43) 

b\ = R\C\ + T 2 ib 2 , (2.44) 

which can be manipulated such that 

ci = ~ ( 2 - 45 ) 

1 12 1 12 


bi = 


Ri 

T] 2 


c 2 + 


TuT 2 i — R\R 2 

7,2 


(2.46) 
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These equations lead to the following matrix expression 


Cl 

bi 




(2.47) 


This relationship describes the forward and reverse traveling waves immediately to 
the left and right of the interface. With the region assumed to be linear, homogeneous 
and isotropic, the waves c 2 and b 2 can be simply related to the corresponding waves 
C 2 and b 2 located a distance l from the interface as 


c 2 = c' 2 e ll 


b’ = b 2 e~ l1 


(2.48) 

(2.49) 


This relationship can further be represented in matrix form as 


e< 1 0 

0 e 7 " z 


c 2 

b- 2 


(2.50) 


The A-parameter relationship between waves ci, b\ and c 2 , b 2 is as follows 


1 ( e< 1 R ie ~ l1 

T i2 \ R ie ii e-'ri 


c 2 

b- 2 


A\\ A12 \ / C 2 

A 2 1 A 22 ) \ b 2 


(2.51) 


Thus, the generalized A-parameter wave matrix can be represented as 


/ An 

Ai 2 

_ 1 [ 

l A 2 i 

a 22 ) 

1 + R 1 


zA Re" 7 * \ 1 / 1 —R 


1 + R \ Re 11 e~A 1 ~ R \ -R 1 


(2.52) 
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This equation may be simplified by letting P = e ll and applying matrix multiplica¬ 
tion to produce 

( A “ ^ 'l = 1 ( 1 - R2pl - P ! ) \ (2 53, 

[a* a 22 ) p(i-ip){r(i-p*) p’-jp ) 1 ' ; 

The Scattering parameters can be determined from the well known relations 


/ S n 

Sl2 > 

1 

A 2 1 AnA 22 — A21A12 j 

[ s 2 1 

*5*22 ; 


v 1 — A 12 / 


Through substitution, you obtain the scattering relations 

Sn=S 22 = R 1 ^~J^ ( 2 . 55 ) 

S 21 =S l2 = ^=P± ( 2 . 56 ) 

With an analytic formulation of the theoretical scattering parameters developed, the 
general process of the Newton-Raphson method can then be invoked to determine 
the complex roots of the nonlinear function, which in this case produces the relative 
permittivity of the measured sample. Note, if the sample is also magnetic, a two- 
dimensional root search can be used to extract both e and //. In this case, Sn and 
5*21 are required. 


2.10 Numerical Position Independent Method 

Simple techniques exist to de-embed the sample position planes from the cali¬ 
brated reference planes, though this requires exact knowledge of the sample position. 
Therefore, a formulation independent of the sample position would lead to greater pre¬ 
cision of the extracted material parameters. This formulation can be accomplished 
using a computational algorithm, such as the Newton-Raphson method, to derive 
complex material parameters from the measured scattering parameters. 


21 



The position independent formulation effectively eliminates any uncertainties 
in the sample position. The Newton-Raphson method can be modified to position- 
independently find the complex parameters of a dielectric and magnetic material 
[8]. In this case the method must be adapted to solve for the complex roots of two 
nonlinear coupled equations, such as 


f(x,y) = \S t n{uj,e,^)S t 2 2 y (uj,e,fi) - S™ P S 22 P \ < tol (2.57) 

g(x, y) = |S^ w (w, e, e, y) - Sf^S^ | < tol (2.58) 

With position uncertainty completely removed, the two functions f(x,y) and g(x,y) 
can be shown in Taylor expansion form 

f(x + Ax,y + Ay) = f(x, y) + x!A — (2.59) 

g(x + Ax,y + Ay) = g(x, y) + + S ^^-Ay. (2.60) 

ox oy 

This can be manipulated into matrix format as 


fx 


9x 



(2.61) 


where the application of the inverse matrix will isolate the new guesses for the two 
coupled nonlinear equations as 


I x i+1 
\ Vi +1 


Xi 


Vi 


1 

fx9y fy9x 


9y fy 

9x fx 



(2.62) 


The numerical partial derivative of f(x) and g(x) can be approximately computed as 

f(x + Ax) - f(x) f(x + Ax) - f{x) 


f x (x) = lim 
Jxy ’ 


Ax 


Ax 


(2.63) 
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g x {x) = lim 

A.x.^0 


g(x + Ax) — g(x) g(x + Ax) — g(x) 


(2.64) 


Ax Ax 

where the partial derivatives with respect to x can be computed in similar format. The 
general process of the Newton-Raphson method then follows where the subsequent 
iterations determine the complex roots of the nonlinear functions, which in this case 
produce the relative permittivity and permeability of the measured sample. 


2.11 Summary 

A discussion of permittivity and permeability was necessary as the extraction 
of these parameters are the main focus of this research. A brief description of the 
Low Frequency Focus Beam System was given as well as a break down of the com¬ 
ponents that allow this system to be a viable way of measuring material constitutive 
parameters. A review of calibration techniques that are in practice now at GE for the 
free space measurement system was provided. This research leads to the conclusion 
that there is not one overall superior calibration technique, they all have positives 
and negatives. Bandwidth limitations and calibration standard errors were two of the 
main limitations of most techniques. There is no way in practice to directly extract 
material constitutive parameters, so various material parameter extraction techniques 
were discussed. These various techniques should help in the analysis of the three short 
calibration technique as information gained from the reflective standards may prove 
to be critical. A discussion of the current time-gating based calibration technique 
used by General Electric was given in an effort to exemplify the need for an improved 
calibration technique. 
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III. Development of the Three Short Calibration Technique 


I n the previous chapter, it was shown that an alternative focused beam calibration 
technique may be beneficial for extracting meaningful uncorrupted data near band 
edges. In this chapter, a three-short calibration technique which does not rely on 
windowing or gating is proposed and developed for the focused beam system. 

3.1 Universal 12-Term Error Model 

An effective method to calibrate the Low Frequency Focus Beam System is ex¬ 
plored here, namely, the three short calibration technique. The three short calibration 
technique involves making 3 initial measurements of highly reflective standards placed 
at different specified positions, along with the general “Thru” measurement in order 
to obtain sufficient data to be able to solve for all the error terms that are encompassed 
in the “12-term error model” [18]. The three short calibration procedure involves a 
12-term error model (6-term model forward, 6-term model reverse). By solving for all 
the error terms, an accurate and robust measurement of the material properties can 
be performed. This section will focus on the development of the 12-term error model 
and the three short calibration technique used to solve it. 

The three short calibration technique has been successfully applied to network 
analyzers that have been used in conjunction with stripline systems in the past [9]. 
The motivation for developing theory to apply this calibration technique to the Low 
Frequency Focus Beam System was discussed previously. 

The focus beam system that is available at the Air Force Institute of Technology 
is designed to record data in the ranges between 4-18Ghz. The focus beam system 
used for the experimentation in this thesis available at GE allows for measurements 
to be taken in the ranges between 0.5 and 6GHz, with interest in the LIHF and 
L band frequency range [14], With proper calibration, this allows the engineer to 
accurately predict how a material will interact when illuminated by an electromagnetic 
wave produced by a VHF or UHF radar. VHF and UHF radars are generally used 
for very long-range surveillance, because of the ability to obtain very large power- 
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Forward Model 



Reverse Model 

Port 1 Port 2 



Figure 3.1: This figure represents the forward and reverse signal flow diagram of 

the twelve term error model. 

aperture products at the lower frequencies and their favorable clutter and propagation 
characteristics [14]. With the presence of these threats, it is important to know how 
a material will react in the presence of that field in those low frequencies. 

The relationship between the measured S-parameters at the NWA detector and 
the sample-plane S-parameters defined by the setup of the system is the foundation of 
the 12-term error model. A signal flow diagram is presented in Figure 3.1 in order to 
represent the error terms of the twelve term error model which are defined as follows 

• Ed = forward directivity 

• Es = forward source match 
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• E rt = forward reflection tracking 

• E l = forward load match 

• Ett = forward transmission tracking 

• E x = forward isolation 

• E' d = reverse directivity 

• E' s = reverse source match 

• E' rt = reverse reflection tracking 

• E'j = reverse load match 

• E' TT = reverse transmission tracking 

• E' x = reverse isolation 

The equations to derive the actual sample S-Parameters at the reference planes 
(Sim, S'i 2 A, S 21 A, and S 22 A) as depicted in Figure 3.1 from the measured detector 
S-parameters (Sn ms , S V2ms , S' 2 i ms , and S' 2 2 ms) described by [22] and depicted in Fig¬ 
ure 3.2 are listed below. Notice that each actual S-parameter is a function of all 
four measured S-parameters. The network analyzer must make a forward and reverse 
sweep to update any one S-parameter. 


SuA = 


Slims Ed 

Ert 


1 + S22 ™r E ' D E'< 


Er 


S21^s-E x 


Sl2ms-gS 


1 + E s ) (1 + S22 ^- e ' d e', 


E'rEr 


S2i m s-E x 


Si 2 ms-E'} 


(3.1) 


S 2 1 A — 


S 2 l,ns~E x 


1 + 


S22ms—E'z 


(■ E’s - E l 


i + 

Ejrt 


1 + 


S22m~E’ D 
E' RT 


E'< 


E'lEl 


S 2 lms Ex 


Sl2ms-E' } 


(3.2) 
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S-IOA — 


Sl2ms-E' a 


1 + 


Sllms-E D 

Ert 


(E s - E' l ) 


1 + 


Slims -E d 


E f 


1 + 


S22ms-E' 1 


E'< 


E'lEl 


S21 m s Ex 


Sl2rns-E>, 


(3.3) 


S22A ~ 


0 + SM W^E' s ) - e l ( 


Sl2ms-g', 


1 + 


S22ms-E' 1 


E'c) -E'r.Er 


Silms-Ez 


Sl2m 8 -E'j 


(3.4) 


Solving for the unknown error terms requires twelve measurements [9]. Each 
measurement represents a separate, unique equation; thus, 12 equations are required 
to solve for 12 unknowns. Two equations come from measuring a matched load on 
each cable leading to the DUT: 


SJjf = 0, and 


Four equations come from an empty focus beam system measurement: 


(3.5) 

(3.6) 


S e u = 0, (3.7) 

S 2 ‘i = 1, (3.8) 

S 22 = 0, and (3.9) 

Si = I- (3.10) 


The remaining six equations are from inserting a short in the sample holder. By 
varying the placement of the short, six unique equations (i.e., cal standard definitions) 
can be formed: 

S[f = - e ~ j2k ^- d \ (3.11) 
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Figure 3.2: A depiction of the phase shift caused by the sample width 


nshl —j2kQ(d—w) 

°22 — e ’ 

(3.12) 

nsh2 i 

On — 1, 

(3.13) 

nsh2 _ -j2k 0 (-w) 

°22 — e ? 

(3.14) 

S°f = -e~ j2k ° d , and 

(3.15) 

jsh3 _ -j2k 0 (-[d+w]) 

'*22 — e 

(3.16) 


After the calibration is performed, post processing still has to be implemented 
to take into account the phase shift caused by the sample width d. From Figure 3.2 
the following actual sample S-Parameters are derived. 


ns _ nmeas j2kol a 

°11 — e 


qs _ nmeas j2k$(l a -\-lb~\ — d) 

°21 — °21 e 5 


























S s 22 = S™ eas e j2ko{ - lb ~ d \ and 

ns _ nmeas j2ko(lb—l-\-l a ) 

°12 ~ °12 e 

With the assumption that l„ and lb are both equal to 0, meaning the calibration 
measurement reference plane now is positioned at z — 0, the following relations exist. 


ns nms 

°11 ~ °11 

(3.17) 

ns nms —jkod 

°21 ~ °21 e 

(3.18) 

ns nms —jkod 

°12 — °12 e 

(3.19) 

qs nms —j2kod 

^22 ~ ^22 e 

(3.20) 


With the S-parameters of the actual sample written in terms of the measured 
S-parameters, the material properties can be calculated using one of the various an¬ 
alytical or numerical techniques. 

3.2 Derivation of Short Placements 

The previous section provided a review of the 12-term error model that describes 
the systematic errors that are inherent in network analyzers. The three short cali¬ 
bration technique allows for six of the twelve error terms to be corrected by unique 
placement and measurement of a highly reflective standard during the calibration 
process. This section will focus in detail on how these unique locations are derived. 

The three short calibration technique is a band limited calibration technique. 
The reference planes that are used for direct placement of the highly reflective stan¬ 
dards used to record phase changes are dependent upon the bandwidth that will be 
used in the experiment. When placing the shorts in the system, spacing is extremely 
important. The purpose of measuring a short in three different positions is to gain 
the six unique equations (three from forward measured and three from reverse mea¬ 
surements) required to solve the 12-term error function. If placed improperly, 
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Short 1 


Short 2 


Short 3 



d 

Figure 3.3: A depiction of the translation of the short standards 

and S(i 3 (defined by the Port 1 measured reflection coefficients of short 1 and short 
3 respectively), for example, could be equal (i.e. degenerate). This would result in 
one fewer equations than unknowns. With short one placed at —d/2 and later short 
three placed at d/2 the distance between the two front faces is d as shown in Figure 
3.3. If d is equal to | (A being the wavelength), the effective phase delay seen at each 
short will be the same. In order to correct for this phase ambiguity, one may impose 
an 80% safety margin on that phase limit to ensure uniqueness. Thus, the goal is to 
keep the distance d to less than |A. 

The focus of this thesis is to be able to extract accurate and dependable low 
frequency material constitutive parameter data. As mentioned previously, the three 
short calibration technique is band limited so an initial determination of the band¬ 
width to be used in the experiment is of importance. A bandwidth of 500MHz to 2GHz 
is the chosen focus for all the experimentation. With an upper frequency established, 
derivation of appropriate spacing may be initiated. 
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The first step in this process is to calculate the |A distance for maximum short 
spacing. This derived distance will nullify the possibility of one short being a degen¬ 
erate of another. The maximum frequency for this experiment is established to be 
2GHz. By placing this frequency into the equation for wavelength, it is shown that 


A — — — 15 cm 


(3.21) 


Dividing this wavelength by 2 and multiplying by the 80 percent safety factor estab¬ 
lishes the maximum spacing distance d between the shorts which is shown to be 


15 cm 


= 6 cm. 


(3.22) 


By specifying a maximum spacing distance of 6cm, instead of 7.5cm (without multi¬ 
plication by the safety factor), one can theoretically evade the possibility of having a 
degenerate error. 


3.3 Phase Angle of Multiple Shorts 

There is another alternative viewpoint that can be implemented to ensure that 
the distance between the shorts will not introduce non-uniqueness into the results. 
This can be done by examining the phase term of each short and finding its location 
on the unit circle as depicted in Figure 3.4. Upon examining the equation defining 
Sfl ' 1 , one can easily deduce the phase term as 

@shi = 2A) 0 d (3.23) 

where 

k 0 = y. (3.24) 

In order to gain a visual perspective behind the phase term, 2n will be represented as 
360°. The wavelength at 500MHz is 60cm long and the wavelength at 2GHz is 15cm. 
With all of these constants defined, a mathematical check of the difference in phase 
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from one translation to the next may be shown. 


2n 360° 

@ s h i = 2k 0 (3cm) = 2 x —3cm = 2 x - x 3 cm = 144° 

A 15cm 

2v T 360° 

©sh 2 = 2/c 0 (Ocm) = 2 x — x 0 cm = 2 x - x 0 cm = 0° 

A 15cm 

2tj- 360° 

@ s /j 3 = 2k 0 (—3cm) = 2 x —d = 2 x -x (—3cm) = —144° 

A 15cm 

The phases 144°, 0°, and —144° are very well spaced around the unit circle. This gives 
confidence that at 2GHz, the phases of the uniquely placed shorts will not introduce 
error into the results due to possible degeneracy. With the upper limit theoretically 
tested, the same process must now be completed for the lower limit of 500MHz to 
ensure sufficient uniqueness of the short standards. Analyzing the multiple shorts at 
500MHz produces the results 

2tt 360° 

Q s hi = 2k 0 (3cm) = 2 x — d = 2 x - x 3cm = 36° 

’ A 60cm 

2t t 360° 

©s /12 = 2h o (0cm) = 2 x —cl = 2 x - x 3cm = 0° 

° v ; A 60cm 

2tj- 360° 

@ s /) 3 = 2ko(—3cm) = 2 x —d = 2 x -x (—3c)m = —36° 

A 60cm 

Although the difference between phases decreased with decreasing frequency, the ob¬ 
tained phase angle difference of 36° is sufficiently resolved by the NWA to ensure 
uniqueness. 

3-4 Time Delays 

In order to properly calibrated the VNA, the exact placement of the shorts 
have to be recorded in terms of the time delays with respect to the z = 0 reference 
plane, with z = 0 being the position where the sample will be placed. In addition the 
Empty/Thru time offset will be set to 0. 
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Figure 3.4: This figure represents the phase angles derived from the placement of 
the multiple shorts. The black marks represent the phase angles obtained at the high 
end of the bandwidth (2GHz) and the red marks represent the phase angles obtained 
by the low end of the bandwidth (500MHz) 

In referencing the multiple shorts by measurements made from Port 1 (forward 
direction), the time constant offsets are defined as 


where d 
light, c, 


t f ~ 

u sh2 — 


-d 

c 

(3.25) 

-, and 

(3.26) 

c 

d 

(3.27) 

c 

= 0 reference plane (3cm) 

and the speed of 


is a constant measured at 299,792,458 meters per second. 
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Table 3.1: This table represents the time delays that were derived and input as 

constants into the network analyzer. The subscripts / and r representing forward 
and reverse direction respectively. 


Short 

Time Delay (ps) 

Shl f 

-100.0692 

Shl r 

89.47857 

Sh2 f 

0 

Sh2 r 

-10.59066 

Sh3 f 

100.0692 

Sh3 r 

-110.65989 


In referencing the multiple shorts by measurements made from Port 2 (reverse 
direction), the time constant offsets are defined as 


f r 

Lch 1 - 


t r — 


d — w 

c 

(3.28) 

— ,and 
c 

(3.29) 

— (d + w) 

(3.30) 


where w is the width of the reflective standard. This width must be taken into account 
in order for accurate time constants to be defined. 


Table 3.1 displays the derived time constants that will be input into the network 
analyzer. Definition of these time constants solves for six of the 12 equations in the 
universal 12 term error model. As previously discussed, four more equations will come 
from an empty/thru measurement where tdeiay will be set to 0. The last two equations 
come from the Isolation measurement, which is used to calibrate out the cross talk 
terms associated with the 12 term error model. With all these measurements in place, 
the NWA is now ready for material measurements. 


3.5 Imprecise Translations 

A cause of error in the results obtained is may be attributed to improper po¬ 
sitioning of the short reference planes. In order for this technique to be effective, 
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precise placement of the shorts is needed. The translation mechanism used for the 
LFFB is shown in Figure 3.5. This mechanism appears to be accurate down to .001 
inches. The reference plane is translated by rotating this crank until the new posi¬ 
tion is reached. There is not only the obvious chance of mechanical error associated 
with this mechanism, but also the added chance of human error in operating the 
mechanism. Although care was taken in all measurements, there was not a way to 
truly validate precise translation. Figure 3.5 shows how the 3-Short Positions were 
marked on the test bench. These positions were not only attempted to be reached 
by the translation mechanism, but also they were measured with a ruler and marked 
as depicted. The ruler measurements and the translation mechanism did not agree 
in all translations, but accuracy was assumed to be within ± .005 inch. Short 1 was 
measured at the mark farthest to the left, followed by Short 2 in the middle (z = 0), 
and Short 3 farthest to the right. Once all 3 Shorts were measured, the sample holder 
was returned to the position of Short 2 at z = 0 and sample measurements were 
taken. One can easily see that without a very precise way to define the reference 
planes, the multiple translations can easily cause unwanted phase errors, and in turn 
add error to the results. In future works, this can be quantitatively validated by per¬ 
forming an uncertainty analysis, ft is noted that a precision laser alignment system 
be implemented in future works. 

3.6 Summary 

The well known error terms that are associated with network analyzers were 
reviewed and discussed. Characterization of these error terms are sought after by 
the use of calibration procedures. The three short calibration technique allows twelve 
measurements to be performed in order to solve for the twelve unknown error terms. 
Since the three short calibration technique is band limited, a derivation of the max¬ 
imum spacing constraint was established and a distance of ±3cm from the z = 0 
reference plane was computed for the short placements. This distance includes a 20% 
safety factor that ensures non-degeneracy. The phase angles produced by the place- 
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merit of the shorts was investigated. The maximum frequency of the bandwidth (i.e. 
2 GHz) produced a phase angle difference of ±144°, between shorts and a phase angle 
difference of ±36°, at 500MHz, which is well within the phase accuracy of the NWA. 
The appropriate time constants were input into the NWA calibration kit definition 
and a discussion of phase advances and phase delays was provided. In the next chap¬ 
ter, the extracted results of two dielectric materials will be presented based on the 
time-domain and the proposed three-short calibration technique outlined above. 
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Figure 3.5: Translation Mechanism 
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IV. Results 


I n this chapter, the permittivity of two dielectric samples will be extracted based on 
the time domain and the proposed three short calibration technique. All measure¬ 
ments were performed using the Low-Frequency Focus Beam System in the frequency 
regime of 0.5-2GHz. Table 4.1 presents the test matrix used for experimentation in 
this thesis. In each measurement, all four scattering parameters were recorded. It 
should be noted that a time domain uncalibrated metal plate response was recorded 
to give physical insight to the system imperfections. This was proceeded by the cal¬ 
ibration standards,three shorts and a THRU measurement, with their results being 
used to solve the 12-term error model. Lastly, calibrated measurements were taken of 
the two known dielectrics. It should be noted that the isolation measurements were 
omitted under the assumption that cross talk was negligible. 


Table 4.1: This table presents the test matrix that will is used during experimen¬ 
tation. The position z = 0 refers to the desired reference plane where actual sample 

measurem ents will be taken. ___ 

Sample Holder Position Object Width(inches) 

z — 0 (uncalibrated) metal plate 0.125 

BEGIN CALIBRATION 


3cm toward Port A (SHORT 1) 

metal plate 

0.125 

z = 0 (SHORT 2) 

metal plate 

0.125 

3cm toward Port B (SHORT 3) 

metal plate 

0.125 

z = 0 (THRU) 

empty 

0 


APPLY CALIBRATION 



4-1 Time Domain Response of Metal Plate 

In order to do a complete analysis of the Three Short Calibration Technique as 
applied to the LFFB system, multiple measurements were taken. The Erst measure¬ 
ment that will be examined is the time domain uncalibrated response of the metal 
plate used as the reflective standard during the calibration process. Viewing this 
data in the time domain allows for the engineer to physically see the obstructions in 
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0.25 



Figure 4.1: Time domain response of metal plate at sample measurement plane. 

the system that may cause interference, which would lead to ripples in the extracted 
constitutive parameter data. 
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In viewing Figure 4.1, one can easily note the main responses of the system. 
At around 27 nanoseconds, the first response occurs. It is well known that electro¬ 
magnetic waves travel at approximately 1 foot per nanosecond in free space. We 
must also note that this is two way time travel, so a response at 27 nanoseconds, is 
indicative of a reflection at approximately 13.5 nanoseconds one way travel, or 13.5 
feet from the NWA. It must be noted that this is only an approximate location due to 
the fact that electromagnetic waves travel at a different speed when in coaxial cables, 
antennas, and dielectric lenses. Included in this length of travel are cables leading to 
the feed horn, distance through the horn, and the free space propagation until the 
wave strikes the front interface of the collimating lens. Also a second reflection occurs 
in this figure. This reflection is the response of the wave traveling through the front 
face of the lens, striking the back face at the exiting point of the lens, and reflecting 
back through the lens to the NWA. The third major reflection that can be noticed 
from this figure is the response from the metal plate, located at approximately 43 
nanoseconds. One can further conclude that these components are indeed the causes 
of the responses by noting that there is no response coming from the second lens. All 
of the energy is reflected back by the metal plate and never reaches the second lens. 
In trying to solely focus on the response of a device under test, the response as a 
result of the reflection from the lens is a noticeable problem. The air/collimating lens 
match seems questionable and could use further investigation. 

There are also multiple responses occurring in the data after the main response 
of the metal plate. This suggests that there could be a coupling issue with the system 
caused by ringing in the lenses or even multiple bounces between the lenses and the 
actual sample. Coupling issues are normally suggested to be mitigated using a time 
domain gating process [11]. However, a deeper concern is system multipath return 
that appears at the same time as the desired metal plate response. It is critical to 
note that this undesired multipath response cannot be gated out, and can therefore 
lead to gross-errors. 
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The distance between the focusing lens and the sample holder is reported to be 
26.75 inches. With the mutipath effects that are noted in Figure 4.1, one can conclude 
that this may be an issue. At the lowest frequency (500MHz) the wavelength (23.6 
inches) begins to approach the actual distance between these two components. At 
this wavelength, it may become electrically difficult to be able to distinguish between 
the two point scatterers, which could lead to a sample-lens coupling effect. 

4-2 Time Domain Gated Data 

As previously discussed, windowing causes data corruption in the low frequen¬ 
cies. To illustrate the need for a better approach, time domain gated data was taken 
for a fiberglass panel. The frequency range of this data is the full range of the LFFB 
system, 500MHz to 6GHz. A 2 nanosecond gate was used to mitigate any reflections 
outside of the plate reflection response. As can be noted in the upper plot of Figure 
4.2, there is a high degree of variation present in the lower frequency range. The 
cause of this data corruption is commented on in Chapter 2. With the application 
of the three short calibration technique, some of the low frequency data should be 
recoverable. Also there is a roll off effect at the tail end of the data due to the data 
lost while applying the window in the frequency domain. 

For completeness another set of windowed data was taken to show the effects 
of windowing in the extracted constitutive parameter data. In the bottom plot of 
Figure 4.2 a 1.25 nanosecond gate was applied. This set of data seems to be a little 
cleaner than the wider gate of 2 nanoseconds, but the obvious corruptions are still 
apparent. There are oscillations beginning from the low end of 500MHz and seemingly 
not settling until 2.3GHz. The roll off effect at the high frequency is still obvious. 

4-3 Calibrated Constitutive Parameter Values of the Fiberglass Panel 

With an understanding of the causes of system corruption and error, a detailed 
look at the calibrated results of a fiberglass panel may now be analyzed. The sample 
under test is 0.183 inches thick, ft is a dielectric with an expected dielectric constant 
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Figure 4.2: TE Incidence Permittivity Measurements. The top plot represents a an 
application of a 2 nanosecond gate. The bottom plot represents an application of a 
1.25 nanosecond gate. 

as noted from the time domain gated data of approximately 4.8. Plots of multiple 
material parameter extraction techniques are presented as a means of comparison. 

The different approaches in calculating the material constitutive parameters will 
allow for better speculation of error within the data. Located on the plots are the 
derived data. These plots are created by use of the 1-Dimensional Newton-Raphson 
Root Search Algorithm, under the assumption that the material is non-magnetic. 
Theoretically, the S 21 and S 12 data should give better results because they are based 
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on information gained from material transmission data. This is due mainly to the fact 
that the precise location of calibration planes is extremely crucial. After translating 
the reflective standard during the calibration process, one can induce error into the 
results by not returning the sample holder to the initial z = 0 calibration plane. 
Investigation of this hypothesis is possible through analysis of the plots presented in 
Figure 4.3. 

The next two plots depicted in Figure 4.3 are data extracted using the Nicolson 
Ross Weir Technique discussed in Chapter 2. One plot is generated by using the 
forward parameters £21 and An, and the other by using the reverse parameters Si 2 
and A 22 . This technique uses both reflection and transmission coefficients in order 
to extract the material constitutive parameter data. This technique may prove to be 
not as effective with this type of system because of the use of the reflection coefficient 
data. 

Also, there are plots present to represent extraction by use of the 2-Dimensional 
Newton-Raphson Root Search Algorithm. This algorithm assumes no known knowl¬ 
edge of whether the material is a dielectric or not. They are added for completeness 
and should be noted that it also uses reflection coefficient data. 

A numerical position independent complex material parameter extraction tech¬ 
nique plot is also presented. This particular numerical method is an adaptation of 
the complex Newton-Raphson root search. The position independent formulation ef¬ 
fectively eliminates any uncertainties in the sample position, which may prove to be 
the best parameter extraction technique. 

In Figure 4.3, there are several details to be noted. The first detail is the 
corrupted data present in the 0.5GHz to 1GHz region. An aim of the three short 
calibration technique applied to this particular system was to be able to recover as 
much low end bandwidth data as possible. It is obvious that there are certain system 
imperfections in the low frequencies that were unable to be calibrated out with the 
use of this technique. 
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Frequency (GHz) 


Figure 4.3: Real and Imaginary Extracted Components of the Fiberglass Panel 

Relative Permittivity 

The cause of this can be related back to Figure 4.1. The initial reflections 
caused by the collimating and focusing lenses at these frequencies are calibrated out, 
but a coupling effect caused by the electromagnetic wave bouncing around within the 
lens was unable to be mitigated. The beam waist at these lower frequencies is also 
extremely wide, possibly causing other obstructions to be illuminated other than just 
the pure sample. The feed horns at these low frequencies also begin to perform rather 
poorly due to wide bandwidth. The data set that seems to have the lowest amplitude 
of variation is the position independent method. 


44 

















As the frequency increases, the results become better. All of the different pa¬ 
rameter extraction techniques seem to merge together and they converge towards the 
actual sample permittivity value. The noise that is visually still present in these fre¬ 
quencies could be due to a couple of factors. The first and most serious factor is that 
a high degree of accuracy when translating the three shorts is needed. The beam 
waist and poor performance of the horn can also be possible sources of error. Phase 
ambiguity may play a role in some of the lingering oscillations still present at the 
upper end of the bandwidth, but it is concluded that it is not a dominate source of 
error due to the fact that the results improve with increasing frequency. 

The next plot to analyze is the imaginary component of the fiberglass panel. 
This panel is a dielectric and is assumed to be lossless. As depicted in the second plot 
of Figure 4.3, the imaginary component definitely converges to the expected value. 
Once again the position independent method seems to be the most accurate material 
parameter extraction technique for this system. 

For completeness, the permeability of the dielectric panel is plotted. These 
results once again agree with our expected results. Once again as frequency increases, 
the parameter extraction techniques tend to converge to the same value. It is worth 
mentioning how the position independent method seems to be an averaging value 
of forward and reverse parameters. This phase shift between forward and reverse 
parameters is extremely evident in the lower frequencies. This can be caused by the 
sample holder not being returned with precision to the z — 0 plane. Nevertheless, 
the position independent extracts a value of 1 which is expected. The imaginary 
component of the permeability is 0, as anticipated 

4-4 Calibrated Constitutive Parameter Values of the Plexiglas Panel 

This section involves a similar analysis of a Plexiglas panel measured via the 
LFFB system. This Plexiglas panel is another known lossless dielectric, therefore the 
imaginary components are expected to be 0 and the theoretical relative permeability 
value should be 1. 
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Figure 4.4: Real and Imaginary Components of Fiberglass Permeability 

Once again, it can be concluded for the sample that the coupling effects of 
multiple bounces within the lens and lens/sample interaction are two notable causes 
of significant error. The position independent method proves to be the most accurate 
of all the material parameter extraction techniques. As noted in the analysis of the 
fiberglass panel, all of the material parameter extraction techniques converge towards 
a common value as the frequency increases due to better system performance and less 
corrupted data in the reflection coefficients. 

It is also worth mentioning that in the lower end of the bandwidth, the results 
obtained as compared to the fiberglass panel are less noisy. The Plexiglas sample 
theoretically has a lower dielectric constant than the fiberglass panel [2], This means 
that the Plexiglas sample is more transmissive, therefore it should have less of a 
coupling effect with the focusing lens due to a smaller amount of energy being reflected 
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Figure 4.5: Real and Imaginary Components of Plexiglass Permittivity 

back. The results obtained confirm this, and suggest that a further investigation on 
the lens/sample interaction would be of interest. 

4-5 Phase Ambiguity 

The theory developed in Chapter Three suggests that there is a maximum spac¬ 
ing constraint that should help to eliminate degenerate errors in the calibration pro¬ 
cess. At the low end of the bandwidth, the smallest change in phase occurs, with a 
change of ±36°. This is theoretically an acceptable value and in reviewing the results 
obtained, it can be confirmed that phase ambiguity is not as big of a concern as com¬ 
pared to the inherent low frequency coupling effects of the system. If phase ambiguity 
was an issue, one would expect for the results to become more unstable as frequency 
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increases. This is clue to the fact that with increasing frequency, wavelength becomes 
shorter and more precision is needed with shorter wavelengths than with longer ones. 
The results presented in Figures 4.3 and 4.5 show that phase ambiguity is not the 
main concern. The extracted data is better with increasing frequency, suggesting 
other causes of error. 

4-6 Sample-Lens Coupling 

With the results obtained not being as smooth as theoretically predicted, it 
is advantageous to explore the possibilities of the sources of error present. Taking 
a closer look at this data before it is processed through the extraction algorithms 
will give insight if the calibration was successful or not. In order to do this, the 
individual S-parameter components were separated. Similar to the process of time 
domain gating, an inverse Fourier transform was applied to the the calibrated S\ \ 
and 5*21. The Matlab fftshift command was used to help visualize the transform 
with the 0 frequency component in the middle of the spectrum. In regards to post 
calibration, one would expect the multiple reflections caused by the lens in Figure 
4.1 to be suppressed and the only response to be left is that of the actual sample. 
In reviewing Figure 4.6, it is evident that the multiple reflections noted in Figure 
4.1 are suppressed. The three short calibration technique does an excellent job of 
mitigating all of the corruption that is present up until the point where the actual 
electromagnetic wave is incident upon the sample. However, it is clear that there are 
oscillations present directly after the main response of the actual sample that were 
not able to be calibrated out. 

It can be concluded that these oscillations are the effects of coupling between the 
lens and the actual sample. At low frequencies where the wavelength is long, this can 
easily become a problem and there has yet to be a study to show how to counteract 
this interaction between the lens and the sample in free space systems without the 
use of time domain gating. 
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Figure 4.6: The top graph represents the Su response of the fiberglass sample post 
calibration. The bottom graph represents the S 21 response of the fiberglass sample 
post calibration 


In order to further confirm the hypothesis of coupling within the components of 
the system, an empty (air) post calibration measurement was taken. Theoretically, 
after calibration one would expect an Su response of zero if there is not a device/sam¬ 
ple under test in the sample holder. This is due to the fact that the calibration process 
aims to mitigate all extraneous sources of clutter, effectively isolating the backscat- 
ter from only the device/sample being measured. In the top plot of Figure 4.8, it 
is clearly seen that there is clutter in the lower end of the bandwidth that was not 
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Figure 4.7: The top plot represents the calibrated time-domain Sn response of an 
empty “air” measurement. The bottom graph represents the calibrated time-domain 
5*2, response of an empty “air” measurement. 

accounted for during the calibration process. This clutter directly manifests itself 
within the sample data of interest. As frequency increases, the apparent corruption 
in the low frequencies becomes less apparent. This suggests that there are scatter¬ 
ing phenomenons due to low frequency coupling between the metal plate used for 
calibration and the focusing lens that were not calibrated out. 

Furthermore, one would expect an S 21 response of one under the assumption 
that no energy is lost as the electromagnetic wave propagates through the system from 
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the transmitting antenna to the receiving antenna. In the bottom plot of Figure 4.8, 
the low frequency coupling phenomenon is once again shown. In correlation to the 
S 11 response, this clutter phenomenon becomes less apparent as frequency increases. 
It is important to take a look at the phase shift caused by this low frequency coupling 
effect. In the lower end of the bandwidth as depicted in Figure 4.9, a peak phase 
shift of ±3° is realized. This phase shift realized by an empty “air” measurement is 
an indicator that the sample data in this region will be corrupted. It is not until you 
get above 1GHz that an average phase shift of less than a degree is realized. Sliding 
window statistics are presented in Tables 4.3 and 4.4. These tables give a better 
insight into the clutter caused by this low frequency coupling effect inherent to the 
LFFB system. 

In Figure 4.7, it is clear that the coupling effect caused by the lenses is not easily 
gated out. In both the S 21 and Sn empty responses, there is a ringing effect that is 
not calibrated out that will cause the extracted sample constitutive parameter data 
to be corrupted. 

4-7 Sample Size 

The size of the sample being measured is an important factor. The phase 
difference produced by the sample is a function of its width and the frequency at 
which it is being tested. The resolution of the NWA becomes better as this phase 
difference increases. The phase angle as a function of frequency and sample width 
for both samples are plotted in Figure 4.10 and are derived from Equation (2.28). It 
can be suggested that the phase angles in the lower region of the bandwidth may not 
be enough for the NWA to completely resolve the signal from the noise floor. This 
phase problem theoretically diminishes with increasing frequency as the wavelengths 
become shorter and a higher phase difference is achieved. System coupling is still 
concluded to be the major source of error, but it is recommended that larger samples 
are used in future experimentation. 
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Table 4.2: This table represents the extracted relative permittivity values of the 

fiberglass and plexiglas samples. Average and standard deviation values over 2 band- 
widths were calculated using the position independent method. 



Fiberglass 

Plexiglas 

average e r (0.75-2GHz) 

5.0959 

2.7754 

standard deviation (0.75-2GHz) 

0.6083 

0.3803 

average e r (1 GHz-2 GHz) 

4.9197 

2.6647 

standard deviation (lGHz-2GHz) 

0.3373 

0.2376 


Table 4.3: This table presents the phase statistics of a calibrated “air” measure¬ 

ment. The bandwidth of the experiment was split into 8 bins in order to view the 
average and variation over specific frequencies. 



S21 

S12 j 

Frequency 

mean 

std 

mean 

std 

500-687.5MHz 

-0.1522 

1.2723 

-0.4527 

1.1572 

687.5-875MHz 

-0.3770 

0.6430 

-0.2977 

0.5422 

875-1062.5MHz 

-0.1237 

0.4912 

0.0822 

0.5053 

1.0625-1.25GHz 

-0.0823 

0.3440 

-0.0935 

0.3947 

1.25-1.4375GHz 

-0.0527 

0.2114 

-0.0692 

0.2750 

1.4375-1.625GHz 

-0.0384 

0.2330 

-0.0530 

0.2572 

1.625-1.8125GHz 

0.0190 

0.1601 

-0.0203 

0.1724 

1.8125-2GHz 

-0.0376 

0.1280 

-0.0479 

0.1408 


Table 4.2 provides the extracted permittivity values of the fiberglass and plex¬ 
iglas samples. Average values over two distinct bandwidths are provided to serve as 
a comparison to the extracted time domain gating method presented in Figure 4.2. 
Excluding the lower end of the bandwidth from 500MHz to 749MHz, results obtained 
agree with the theoretical values of the two dielectrics. Standard deviation values 
were calculated to represent the variation in data due to the low frequency coupling 
effects. 
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Table 4.4: This table presents the magnitude statistics of a calibrated “air” mea¬ 
surement. The bandwidth of the experiment was split into 8 bins in order to view the 
average and variation over specific frequencies. 



S2! 

S!2 

Frequency 

mean 

std 

mean 

std 

500-687.5MHz 

0.9931 

0.0210 

0.9948 

0.0192 

687.5-875MHz 

1.0035 

0.0109 

0.9995 

0.0087 

875-1062.5MHz 

1.0056 

0.0081 

1.0073 

0.0067 

1.0625-1.25GHz 

0.9989 

0.0060 

0.9972 

0.0064 

1.25-1.4375GHz 

1.0004 

0.0037 

1.0023 

0.0044 

1.4375-1.625GHz 

0.9993 

0.0041 

0.9992 

0.0045 

1.625-1.8125GHz 

0.9998 

0.0028 

0.9987 

0.0027 

1.8125-2GHz 

1.0004 

0.0023 

1.0009 

0.0024 
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Figure 4.8: The top plot represents a calibrated Sn response of an empty “air” 

measurement. The bottom plot represents a calibrated £21 response of an empty 
“air” measurement. 
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Figure 4.9: This plot represents the phase shift caused by the low frequency coupling 
effect. 
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Figure 4.10: The top plot represents the phase angle difference realized by the 0.183 
inch fiberglass sample. The bottom plot represents the phase angle difference realized 
by the 0.175 inch plexiglas sample. 
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V. Conclusions 


T he three short calibration technique was examined on a free space measurement 
system, the Low Frequency Focus Beam System. Although theory suggests 
that this technique is applicable, the results obtained suggest low frequency cou¬ 
pling between the multiple components of the system cause errors that are not easily 
calibrated out. Different numerical extraction techniques were used as a basis of com¬ 
parison. The extracted constitutive parameters for both the Fiberglass sample and 
the Plexiglas sample, excluding the corrupted lower end of the bandwidth, on average 
were in agreement with the time domain gated data obtained by General Electric. 

The main source of error within the results obtained is concluded to be caused 
by coupling effects within the system. By evaluating a post calibration empty (air) 
measurement, it is experimentally shown that there is a coupling effect associated 
with the dielectric lenses used to collimate and focus the electromagnetic wave onto 
the device under test. It is further shown by evaluation of the Sn and S 22 responses 
of the dielectric samples, that there is also interaction between the sample and the 
lenses. This is validated because of the increased amplitude of oscillations as compared 
to the empty measurement. A position independent material parameter extraction 
algorithm proved to be the most accurate extraction technique as it nullified the phase 
differences between forward and reverse parameters. Although systematic errors were 
encountered, the average experimentally measured material parameter data obtained 
agreed with theoretically documented values. 

5.1 Future Work 

In order to utilize this technique on the Low Frequency Focus Beam System, 
a more precise way of translating the shorts in defining calibration planes is needed. 
Theory suggests that with precise definition of these calibration planes, accurate con¬ 
stitutive parameter data may be extracted. If the bandwidth of interest is 500MHz to 
2GHz, separating the experiment into two bandwidth limits may prove to be useful. 
For example 500MHz-1.25GHZ and 1.25GHz-2GHz. 
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5.1.1 Quiet Zone and Reflective Standard. A study can be conducted on the 
quiet zone of this system. This system is designed to where the energy is focused at 
the focal length of the focusing lens. By translating the short the maximum allowable 
distance, the energy distribution may become corrupted. Also more care can be taken 
when choosing the reflective standard. Literature states that less energy is reflected 
at lower frequencies. Polishing the standard would increase the reflected signal which 
would help separate the received signal from the noise floor of the system. 

5.1.2 Sensitivity Analysis. A quantitative study can be performed on how 
much error is induced into the system through improper calibration plane specifica¬ 
tion. The equations of the universal 12-term error model are well established. Altering 
each error term can be invoked in order to see which error terms cause the most error 
uncertainty within the results. With this information known, a more robust focus can 
be done on controlling the error associated with those terms. 

5.1.3 FDTD Modeling. In order to truly realize the source of the multiple 
reflections, a more extensive time domain modeling of the system should be investi¬ 
gated using the Finite Difference Time Domain technique. 



Appendix A. MatLab Code 


T he Mat lab® code used to produce the plots in this thesis is presented. Multi¬ 
ple material parameter extraction techniques including the NRW and Newton- 
Raphson (ID, 2D, and Position Independent) were discussed and are placed in this 
Appendix for the purposes of reproducing results obtained in this thesis. 


Listing A.l: This hie is the base hie which contains functions to the various material 
parameter extraction techniques 

(appendixl/WGunnrhesf s.m) 

1 

function [] = WGunn_Thesis() 


clear all; close all; clc 

6 length_Sample_plexi =.00439412; 

lengtli_Sample_fg = . 004648199999995 ; 

K_csq =0; 

% Parameter guesses 
11 mu_plexi = 1 - O.OOlj; 

epsilon_plexi = 2.8 - O.OOlj; 

mu_f g = 1 - . 001 j ; 

epsilon_fg = 4.86 - O.OOlj; 

16 

7, Determine the S-parameters for Plexiglass and Fiberglass 
[material_Measurements] = getLabData(); 

Freq = material.Measurements(: , 1) ; 

S21_l = material.Measureients(: ,2) ; 

21 Sll_l = material.Measureients(: ,3) ; 

S12_l = material.Measurements(: ,4) ; 

S22_l = material.Measurements(:,5); 

[material_Measurements] = getLabData(); 

26 Freq = material.Measurements(:,1); 

S21_2 = material.Measureients(:,2); 

Sll_2 = material.Measureients(: ,3) ; 

S12_2 = material.Measureients(:,4); 

S22_2 = material_Measurements(:,5); 

31 


C = 2.997925e8; 

36 mu_0 = pi*4e-7; 

epsilon_0 = 1/(C~2*mu_0); 

K_osq = (2*pi .*Freq). ~ 2 .*epsilon_0.*mu_0; 

gamma_o = sqrt(K_csq - K_osq); 
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41 K_o = sqrt (K_osq); 

S22_1 = S22_1 .* exp (-1i*2 . *K_o*length_Sample_plexi) ; 
S21_1 = S21_1 .*exp(-li.*K_o*length_Sample_plexi) ; 
S12_1 = S12_1 .*exp(-li.*K_o*length_Sample_plexi) ; 

46 

S22_2=S22_2.*exp(-li*2.*K_o*length_Sample_fg); 
S21_2 = S21_2.* exp ( — 1i.*K_o*length_Sample_fg) ; 

S12_2 = S12_2 . * exp ( — 1 i . *K_o*length_Sample_fg) ; 
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•/. Use the NRW method to compute the Plexiglass and Ram material 
parameters 

[mu_rel_forward_1, epsilon_rel_forward_1] = computeNRW(SI1_1, ... 

S21_l , . . . 

gamma.o, length_Sample_plexi, K_csq, K_osq); 

56 [mu_rel_reverse_1 , epsilon_rel_reverse_1] = computeNRW(S22_1, .. 

S12_l , . . . 

gamma.o, length_Sample_plexi, K_csq, K_osq); 

[mu_rel_forward_2, epsilon_rel_forward_2] = computeNRW(SI1_2, ... 

S21_2 , . . . 

gamma.o , length_Sample_fg , K_csq , K_osq) ; 

61 [mu_rel_reverse_2 , epsilon_rel_reverse_2] = computeNRW(S22_2 , ... 

S12_2 , . . . 

gamma.o , length_Sample_fg , K_csq , K_osq) ; 


"/, Use ID and 2D Newton-Raphson complex root search method 
66 [epsilon_S11_plexi] = newtonRootSearchSl1(Freq,length_Sample_plexi 

Sll_l, K_csq, mu_plexi, epsilon_plexi); 

[epsilon_S22_plexi] = newtonRootSearchSll(Freq,length_Sample_plexi 

S22_l, K_csq, mu_plexi, epsilon_plexi); 

[epsilon_S21_plexi] = newtonRootSearchS21(Freq,length_Sample_plexi 

71 S21_l, K_csq, mu_plexi, epsilon_plexi); 

[epsilon_S12_plexi] = newtonRootSearchS21(Freq,length_Sample_plexi 

S12_l, K_csq, mu_plexi , epsilon_plexi) ; 

[epsilon_Sll_fg] = newtonRootSearchSll(Freq,length_Sample_fg , . . . 

76 Sll_2, K_csq, mu_fg, epsilon_fg); 

[epsilon_S22_fg] = newtonRootSearchSll(Freq,length_Sample_fg , . . . 

S22_2 , K_csq, mu_fg, epsilon_fg); 

[epsilon_S21_fg] = newtonRootSearchS21(Freq,length_Sample_fg , . . . 
S21_2, K_csq, mu_fg, epsilon_fg); 

81 [epsilon_S12_fg] = newtonRootSearchS21(Freq,length_Sample_fg , . . . 
S12_2, K_csq, mu_fg, epsilon_fg); 
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[epsilon_f_plexi1 ,mu_f_plexi1] = newtonRootSearch2D(Freq , . . . 

length_Sample_plexi , Sll_l, S21_l, mu_plexi , epsilon_plexi ,.. . 
K_csq) ; 

86 [epsilon_r_plexi1 ,mu_r_plexi1] = newtonRootSearch2D(Freq , . . . 

length_Sample_plexi , S22_l, S12_l, mu_plexi , epsilon_plexi ,.. . 
K_csq) ; 

[ epsilon_f_fgl ,mu_f_fgl] = newtonRootSearch2D(Freq, ... 

length_Sample_fg , Sll_2, S21_2, mu_fg , epsilon_fg,K_csq); 
[epsilon_r_fgl,mu_r_fgl] = newtonRootSearch2D(Freq, ... 

91 length_Sample_fg , S22_2, S12_2, mu_fg , epsilon_fg,K_csq); 

% Use position independent 2D Newton-Raphson method 
[epsilon_plexi2, mu_plexi2] = poslndependent(Freq, ... 

length_Sample_plexi , . . . 

Sll_l, S12_l, S21_l, S22_l, mu_plexi, epsilon_plexi,K_csq); 

96 [epsilon_fg2 , mu_fg2] = poslndependent(Freq, length_Sample_fg , . . . 
S11_2 , S12_2 , S21_2 , S22_2 , mu_fg, epsilon.fg,K_csq) ; 
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/ Display all the plexiglass parameters 
scrsz = get (0,’ScreenSize’); 

hi = figure (’Name’,’Real Component Plexiglass Permittivity ’ , 1 .. . 


Position 1 , . . . 

[1 scrsz(4)/2 scrsz(3)/2 scrsz(4)/2]) ; 


106 plot(Freq, real (epsilon_S11_plexi) , 

, [1,0,0]) ; hold on 

’LineWidth’ 

,1.5, 

’color 

plot(Freq, real (epsilon_S12_plexi) , 

, [ .9,0.1 , 0.2]) ; hold on 

’LineWidth’ 

,1.5, 

’color 

plot(Freq, real (epsilon_S21_plexi) , 

, [ .8,0.2,0.4]) ; hold on 

’LineWidth’ 

,1.5, 

’color 

plot(Freq, real (epsilon_S22_plexi) , 

, [ .7,0.3,0.6]) ; hold on 

’LineWidth’ 

,1.5, 

’color 


plot(Freq, real (epsilon_rel_forward.l),’LineWidth’, 4 ,'color’... 
,[0,0,1]); hold on 

111 plot(Freq, real (epsilon_rel_reverse_l),’LineWidth’,4,’color’... 

, [. 1, .1, .8]) ; hold on 

plot(Freq, real(epsilon_f_plexil) , ’LineWidth’ ,1.5, ’color’ , [0,1,0]) 
; hold on 

plot(Freq, real(epsilon_r_plexil) , ’LineWidth’ ,1.5, ’color’... 

,[.1,.8,.1]); hold on 

plot(Freq, real(epsilon_plexi2) , ’LineWidth’ ,3, ’color’ , [0,0,0]) 
ylim ( [0 15]) 

116 xlabel (’Frequency (GHz) ’ , ’FontSize’ ,14) ; 
ylabel (’Re(\epsilon)’, ’FontSize’, 14); 

legend(’Sll Derived’,’S12 Derived’ , ’ S21 Derived’ , ’S22 Derived’,’.. 
NRW F ’ , . . . 

’NRW R’,’2D Newton F’,’2D Newton R’,’Pos Ind’,3) 

121 h2 = figure (’Name’,’Imaginary Component Plexiglass Permittivity’,’ 
Position ’ , . . . 
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[1 scrsz(4)/2 scrsz(3)/2 scrsz(4)/2]) ; 


plot(Freq, imag (epsilon_S11_plexi) , 

, [1,0,0]) ; hold on 

’LineWidth’ 

,1.5. 

, ’color ’ . 

plot(Freq, imag (epsilon_S12_plexi) , 

, [ . 9,0.1 , 0.2]) ; hold on 

’LineWidth’ 

,1.5. 

, ’color ’ . 

plot(Freq, imag (epsilon_S21_plexi) , 

, [ .8,0.2,0.4]) ; hold on 

’LineWidth’ 

,1.5, 

, ’color ’ . 

plot(Freq, imag (epsilon_S22_plexi) , 

, [ .7,0.3,0.6]) ; hold on 

’LineWidth’ 

,1.5, 

, ’color ’ . 

plot (Freq, imag(epsilon_rel_forward_ 

. 1) , ’LineWidth’ 

,4,’color 


,[0,0,1]); hold on 

plot(Freq, imag(epsilon_rel_reverse_l),’LineWidth’,4,'color’... 

, [. 1, .1, .8]) ; hold on 

plot(Freq, imag(epsilon_f_plexil) , ’LineWidth’ ,1.5, ’color’ , [0,1,0]) 
; hold on 

plot(Freq, imag(epsilon_r_plexil) , ’LineWidth’ ,1.5, ’color’... 

,[.1,.8,.1]); hold on 

131 plot(Freq, imag(epsilon_plexi2) , ’LineWidth’ ,3, ’color’ , [0,0,0]) 
ylim([-5 5]) 

xlabel (’Frequency (GHz) ’ , ’FontSize’ ,14) ; 
ylabel (’Im(\epsilon )’ , ’FontSize’, 14); 

legend(’Sll Derived’,’S12 Derived ’ , ’ S2 1 Derived ’ , ’ S22 Derived’,’.. 
NRW F ’ , . . . 

136 ’ NRW R ’ , ’ 2D Newton F ’ , ’ 2D Newton R ’ , ’ Pos Ind’,3) 

h3 = figure (’Name’,’Real Component Plexiglass Permeability’,’... 
Position ’ , . . . 

[1 scrsz(4)/2 scrsz(3)/2 scrsz(4)/2]) ; 
plot(Freq, real (mu_rel_forward.l) , ’LineWidth’ ,4, ’color’ , [0,0,1]) ; 
hold on 

141 plot(Freq, real (mu_rel_reverse_l),’LineWidth’ ,4,’color’.. . 

, [. 1, .1, .8]) ; hold on 

plot(Freq, real(mu_f_plexil), ’LineWidth’ ,1.5, ’color’ ,[0,1,0]) ;... 

hold on 

plot(Freq, real(mu_r_plexil) , ’LineWidth’ ,1.5, ’color’ , [.1, .8, .1]) ; 
hold on 

plot(Freq, real (mu_plexi2) , ’LineWidth’ ,3, ’color’ , [0,0,0]) 
ylim([-5 5] ) 

146 xlabel (’Frequency (GHz) ’ , ’FontSize’ ,14) ; 
ylabel (’Re(\mu)’, ’FontSize’, 14); 

legend(’NRW F ’ , ’NRW R’ , ’2D Newton F’ , ’2D Newton R’,’Pos Ind’,3) 

h4 = figure (’Name’,’Imaginary Component Plexiglass Permeability’,’ 
Position ’ , . . . 

151 [1 scrsz(4)/2 scrsz(3)/2 scrsz(4)/2]) ; 

plot(Freq, imag( mu_rel_forward.l) , ’LineWidth’ ,4, ’color’ , [0,0,1]) ; 
hold on 

plot(Freq, imag( mu_rel_reverse _1) , ’LineWidth’ ,4, ’color’ .. . 

, [.1 , .1 , .8]) ; hold on 

plot(Freq, imag(mu_f_plexil),’LineWidth’,1.5,’color’,[0,1,0]); ... 

hold on 
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plot(Freq, imag(mu_r_plexil) , ’LineWidth’ ,1.5, ’color’ , [.1, .8, .1]) ; 
hold on 

156 plot(Freq, imag(mu_plexi2) , ’LineWidth’ ,3, ’color’ , [0,0,0]) 
ylim([-10 10] ) 

xlabel (’Frequency (GHz) ’ , ’FontSize’ ,14) ; 
ylabel (’Im(\mu)’, ’FontSize’, 14); 

legend(’NRW F’,’NRW R’,’2D Newton F’,’2D Newton R’,’Pos Ind’,3) 

161 


7o Display all the RAM parameters 

h5 = figure (’Name’,’Real Component Fiberglass Permittivity’,’... 
Position ’ , . . . 

166 [1 scrsz(4)/2 scrsz(3)/2 scrsz(4)/2]) ; 


plot(Freq, real (epsilon_S11_fg) , 
hold on 

’LineWidth’ 

,1.5, 

’color ’ ,[1,0,0]) 

plot(Freq, real (epsilon_S12_fg) , 
, [ . 9,0.1 , 0.2]) ; hold on 

’LineWidth’ 

,1.5, 

’color ’ . .. 

plot(Freq, real (epsilon_S21_fg) , 
, [ .8,0.2,0.4]) ; hold on 

’LineWidth’ 

,1.5, 

’color ’ . .. 

plot(Freq, real (epsilon_S22_fg) , 

, [ .7,0.3,0.6]) ; hold on 

’LineWidth’ 

,1.5, 

’color ’ . .. 


171 plot(Freq, real (epsilon_rel_forward_2) , ’LineWidth’ ,4, ’color’... 
,[0,0,1]); hold on 

plot(Freq, real (epsilon_rel_reverse_2),’LineWidth’,4,’color’... 

, [. 1, .1, .8]) ; hold on 

plot(Freq, real(epsilon_f_fgl) , ’LineWidth’ ,1.5, ’color’ , [0,1,0]) ; . 

hold on 

plot(Freq, real(epsilon_r_fgl) , ’LineWidth’ ,1.5, ’color’ , [.1 , .8, .1]) 
; hold on 

plot(Freq, real(epsilon_fg2),’LineWidth’,3,’color’,[0,0,0]) 

176 ylim( [0 15] ) 

xlabel (’Frequency (GHz) ’ , ’FontSize’ ,14) ; 
ylabel (’Re(\epsilon)’, ’FontSize’, 14); 

legend(’Sll Derived’,’S12 Derived’,’S21 Derived’,’S22 Derived’,’.. 
NRW F ’ , . . . 

’NRW R’ , ’2D Newton F’ , ’2D Newton R’ , ’ Pos Ind’,3) 

181 


h6 = figure (’Name’,’Imaginary Component Fiberglass Permittivity’,’ 
Position ’ , . . . 

[1 scrsz(4)/2 scrsz(3)/2 scrsz(4)/2]) ; 


plot(Freq, imag (epsilon_S11_fg ) , 
hold on 

’LineWidth’ 

,1.5, 

’color’ , [1 

186 plot(Freq, imag (epsilon_S12_fg) , 

, [ . 9,0.1 , 0.2]) ; hold on 

’LineWidth’ 

,1.5, 

’color ’ . .. 

plot(Freq, imag (epsilon_S21_fg) , 

, [ .8,0.2,0.4]) ; hold on 

’LineWidth’ 

,1.5, 

’color ’ . .. 

plot(Freq, imag (epsilon_S22_fg) , 
, [ .7,0.3,0.6]) ; hold on 

’LineWidth’ 

,1.5, 

’color ’ . .. 


plot(Freq, imag(epsilon_rel_forward_2),’LineWidth’,4,’color’... 
,[0,0,1]); hold on 
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plot(Freq, imag(epsilon_rel_reverse_2),’LineWidth’, 4 ,'color’... 

, [. 1, .1, .8]) ; hold on 

191 plot(Freq, imag(epsilon_f_fgl), ’LineWidth’, 1.5, ’color’, [0,1,0]); . 

hold on 

plot(Freq, imag(epsilon_r_fgl) , ’LineWidth’ ,1.5, ’color’ , [.1 , .8, .1]) 
; hold on 

plot(Freq, imag(epsilon_fg2),’LineWidth’,3,’color’,[0,0,0]) 
ylim([-5 5]) 

xlabel (’Frequency (GHz) ’ , ’FontSize’ ,14) ; 

196 ylabel (’Im(\epsilon )’ , ’FontSize’, 14); 

legend(’Sll Derived’,’S12 Derived ’ , ’ S2 1 Derived ’ , ’ S22 Derived’,’.. 
NRW F ’ , . . . 

’NRW R’,’2D Newton F’,’2D Newton R’,’Pos Ind’,3) 


201 h7 = figure (’Name’,’Real Component Fiberglass Permeability’,’... 
Position ’ , . . . 

[1 scrsz(4)/2 scrsz(3)/2 scrsz(4)/2]) ; 


plot (Freq , 

real (mu_rel_forward_2 ) , 

’LineWidth’ , 4, : 

’ color ’ , [0 

,0,1]) 

hold on 







plot (Freq , 

real (mu_rel_reverse_2 ) , 

’LineWidth’ , 4, : 

’color ’ . .. 


,[.l,.l 

,.8]); hold on 






plot (Freq , 

real (mu_f_fgl) 

, ’LineWidth’ ,1.5, 

’color : 

1 ,[0,1,0]) 

; hold 

206 plot (Freq , 

real (mu_r_fgl) 

, ’LineWidth’ ,1.5, 

’color : 

’ , [.1 , .8, . 

1]) ; . 

hold on 







plot (Freq , 

real (mu_f g2) , ’ 

LineWidth 

’ ,3, ’color ’ , [0 

,0,0]) 



ylim([-5 5] ) 

xlabel (’Frequency (GHz) ’ , ’FontSize’ ,14) ; 
ylabel (’Re(\mu)’, ’FontSize’, 14); 

211 legend (’NRW F’,’NRW R’,’2D Newton F’,’2D Newton R’,’Pos Ind’,3) 

h8 = figure (’Name’,’Imaginary Component Fiberglass Permeability’,’ 
Position ’ , . . . 

[1 scrsz(4)/2 scrsz(3)/2 scrsz(4)/2]) ; 
plot(Freq, imag (mu_rel_forward_2) , ’LineWidth’ ,4, ’color’ , [0,0,1]) ; 
hold on 

216 plot(Freq, imag (mu_rel_reverse_2),’LineWidth’ ,4,’color’.. . 

, [. 1, .1, .8]) ; hold on 

plot(Freq, imag( mu_f _fgl), ’LineWidth’,1.5, ’color’, [0,1,0]); hold . 
on 

plot(Freq, imag (mu_r_fgl), ’LineWidth’,1.5, ’color’ ,[.l,.8,.l]); ... 

hold on 

plot(Freq, imag( mu_fg2 ) , ’LineWidth’ ,3, ’color’ , [0,0,0]) 
ylim([-10 10] ) 

221 xlabel( ’Frequency (GHz) ’ , ’FontSize ’ ,14) ; 
ylabel (’Im(\mu)’, ’FontSize’, 14); 

legend (’NRW F’,’NRW R’,’2D Newton F’,’2D Newton R’.’Pos Ind’,3) 
end 
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Listing A.2: This file allows the data obtained from the NWA to be imported into 

mat lab. 

(appendixl / getLabData.m) 

% Network Anaylzer Data Importer 
7, Lt William Gunn 

7. 

7, [Haterial.Measurements] = getLabData (inputFile) imports data ... 
files from 

5 7. the network analyzer in the materials lab into a single matrix ... 
with each 

7. column representing frequency, Sll, S12 , S21 , S22 . The optional... 

argument 

7. FILENAME is a string that contains the path and name of the data... 
file in 

7. the format DRI VE_LETTER : \ DIRECTORY \ ... \ DIRECTORY\FILENAME . CTI . 

If 

7« FILENAME isn’t specified, the user will be prompted to specify ... 
the file 

10 7. with an open-file dialog box. 

function [Material_Measurements] = getLabData(inputFile) 

7.Allow the user to select the files to import if not specified 
15 if nargin ~= 0 

fid = fopen (inputFile) ; 
elseif (nargin= = 0) I I (fid==-l) 

[impfilenm, impfilepath] = uigetfile ({ 1 *.cti’ , ’CTI Files (*.... 
cti) 1 ; . . . 

’*.txt 1 , ’ Text Files (*.txt) ’;’*.* 1 ,’All Files (*.*)’},’... 
Select Data File for Import’); 

20 fid = fopen ([impfi1epath impfilenm]); 

end 

7. Scan the file and extract pertinent information. Setting the ... 
headerlines 

7. eliminates the entire header such that it begins with the ... 
frequency 

25 C_Data = textscan (f id , ’7«s’, ’ headerlines ’ ) ; 

7. This undoes the strange cell-within-a-cell output that textscan ... 
uses 

C_Data = C_Data{:}; 

30 l Splits the column for S-parameters into a real and imaginary ... 
column 

[C_Data(:,l), C_Data(:,2)] = strtok(C_Data,’,’); 

°l Locate the end of the file heading 
data_Begin = strmatch(’VAR_LIST_BEGIN’,C_Data); 
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•/. Locate the beginning of each data column and define column 
length 

Data_Cells = strmatch( 1 BEGINC_Data) ; 

C_Length = Data_Cells(1) - data_Begin - 2; 

7, Preallocate variables for Freq , Sll, S12 , S21, S22 
var_Meas = length (Data_Cells) + 1; 

Measurements = nan([C_Length, var.Meas]); 


45 

7,Fill measurements array with Frequency, Sll, S12 , S21, S22 

Measurements (:, 1) = str2double(C_Data((data_Begin + l):(C_Length +... 

dat a_Begin))) ; 

parfor i = 1: length (Data_Cel1s) 

50 Measurements(:,i+1) = str2double(C_Data(Data_Cells(i) + 1: ... 

Data_Cells(i) + C_Length , 1))+ j .*str2double(C_Data(.. . 

Data_Cells(i) + 1: Data_Cells(i) + C_Length , 2)) ; 

end 

7o Close the data file 
fclose (fid) ; 

55 

7«Pass extracted information out 
Material.Measurements = Measurements; 

end 


Listing A.3; This file uses the Nicolson-Ross-Weir material parameter extraction 
technique 

(appendixl/computeNRW.m) 

1 7. NRW Method Computation 
°l Lt William Gunn 


function [mu_rel, epsilon_rel] = computeNRW(S11_22, S21_12, 
gamma.o, ... 

length_Sample , K_csq , K_osq) 


7. Calculate values for Q 

Q = (S11_22. ~ 2 - S21_12 . ~ 2 + 1) ./(2 . * SI1_22) ; 

11 

7. Calculate R and determine proper root choice 
R_pos = Q + sqrt(Q.“2 - 1); 

R_neg = Q - sqrt(Q.~2 - 1); 

R = (abs(R_pos)<=l).*R_pos + (abs(R_neg)<=l).*R_neg; 

16 

7. Calculate values for P 
P = S21_12./(1 - R.*S11_22) ; 

/ Calculate the sample relative permittivity and permeability 
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21 mu_rel = (-log(P)./(gamma_o. *length_Sample)). * ( (1 + R) ./(1 - R)); 

epsilon_rel = (K_csq - ( log (P) ./length_Sample) ."2) ./(K_osq.*mu_rel. . . 

) ; 

Listing A.4: This file uses the Newton-Raphson ID material parameter extraction 
technique. This technique assumes that the material is either a dielectric or purely 
magnetic. The reflection coefficients are used to extract material constitutive param¬ 
eters. 

(appendixl/newtonRoot SearchS 11. m) 

7, Complex Newton Root Search 
7, Lt William Gunn 
3 7. 

7. The newtonRootSearchlD () function numerically seeks the complex 
7. permittivity for each S-Parameter. 


8 function [ epsilon_S11] = newtonRootSearchS11(Freq,length_Sample ,.. . 
S11_meas , . . . 

K_csq , mu_rel , epsilon_rel) 

7. Initialize constants 
tol = 10~(-7) ; 

13 eps = 1; 

k = 0; 

delta = le-T; 

18 x = mu_rel; 

y = epsilon.rel ; 

dy = y + delta; 

23 


while ((eps > tol) && (k < 100)) 

7, Compute the theoretical S-parameters 

[Sll_22_thy] = computeS 11 Parameter(Freq,length_Sample , x ,y,.. . 
K_csq) ; 

28 [S11_22_thy_dy] = computeSlIParameter(Freq,length_Sample ,x,dy,. . . 

K_csq) ; 

7. Define the standard f unct ions 
F = Sll_22_thy - Sll_meas; 

33 7» Define the derivative functions 

F_dy = S11 _22_thy_dy - Sll_meas; 

dFy = (F_dy - F)./delta; 

38 7. Calculate next permittivity guess 

yn = y - F./dFy; 
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approximate 


error 


43 


48 


7. Find the lowest relative 
eps = min(abs((yn - y))); 

7. Set the new permittivity 

y = yn; 

°/ 0 Calculate the new deltas 
7, delta = -F./dFy; 


guess 


dy = y + delta; 


k = k + 1; 

53 end 


°/ 0 Pass the relative permittivites out for each S-parameter 
epsilon_S11 = y; 


Listing A.5: This file uses the Newton-Raphson ID material parameter extraction 
technique. This technique assumes that the material is either a dielectric or purely 
magnetic. The transmission coefficients are used to extract material constitutive 
parameters. 

(appendixl/newtonRootSearchS21.m) 

/ Complex Newton Root Search 
% Lt William Gunn 

7 . 

4 7, The newtonRootSearchlD () function numerically seeks the complex 
7. permittivity for each S-Parameter. 


function [ epsilon_S21] = newtonRootSearchS21(Freq,length_Sample ,. . . 

S21_meas , . . . 

9 K_csq, mu_rel, epsilon_rel) 

7. Initialize constants 
tol = 10"(-7) ; 
eps = 1; 

14 k = 0; 


delt 

a 

= le-7; 

x = 

mu 

_rel ; 

y = 

ep 

silon.rel ; 

dy = 

y 

+ delta; 

whil 

e 

(( eps > to 


7, Compute the 


&& (k < 100)) 
theoretical S-parameters 
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[S21_12_thy] = computeS21 Parameter(Freq,length_Sample ,x , y ,.. . 
K_csq) ; 

[S2 1 _12_thy_dy] = computeS21Parameter(Freq,length_Sample , x , dy , . . . 

K_csq) ; 

29 

7, Define the standard funct ions 
F = S21_12_thy - S21_meas; 

°/« Define the derivative functions 
34 F_dy = S21 _12_thy_dy - S21_meas; 

dFy = (F_dy - F)./delta; 

7. Calculate next permittivity guess 
39 yn = y - F./dFy; 

7, Find the lowest relative approximate error 
eps = min(abs((yn - y))); 

44 7, Set the new permittivity guess 

y = yn; 

7. Calculate the new deltas 
'/. delta = -F./dFy; 

49 

dy = y + delta; 
k = k + 1; 

end 

54 

°i Pass the relative permittivites out for each S-parameter 
epsilon_S21 = y; 

Listing A.6: This file uses the Newton-Raphson 2D material parameter extraction 
technique. This technique assumes that the user does not know the material consti¬ 
tutive parameters. The reflection and transmission coefficients are used to solve two 
equations with two unknowns to extract material constitutive parameters. 

(appendixl / newtonRootSearch2D.m) 

"/, Position Independent Analysis 
% Lt William Gunn 

7, 2-Dimensional Newton Root Search Method 

4 


function [epsilon, mu] = newtonRootSearch2D(Freq,length_Sample,... 
Sll_meas, S21_meas, x, y, K_csq) 


% Initialize constants 
tol = 10'(-7) ; 
eps = 1; 
k = 0; 
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14 delta = 10"(-7) ; 


dx = x + delta; 
dy = y + delta; 


19 

while ((eps > tol) && (k < 100)) 

1 Get S-paraneters for standard and derivative functions 
[SI l_22_thy] = computeS11 Parameter(Freq,length_Sample ,x , y . 
K_csq) ; 

[S21_12_thy] = computeS21 Parameter(Freq,length_Sample ,x , y. 
K_csq) ; 

24 [ S11 _22_thy_dy] = coiputeS1IParameter(Freq,length_Sample ,x,dy 

K_csq) ; 

[S21_12_thy_dy] = computeS2IParameter(Freq,length_Sample,x,dy 

K_csq) ; 

[ S11_22_thy_dx] = computeS1IParameter(Freq,length_Sample,dx,y 

K_csq) ; 

[S21_12_thy_dx] = computeS2IParameter(Freq,length_Sample,dx,y 

K_csq) ; 


29 '/. Define the standard functions 

F = Sll_22_thy - Sll_meas; 

G = S21_12_thy - S21_meas; 


'/. Define the derivative functions 
34 F_dy = S11 _22_thy_dy - Sll_meas; 

G_dy = S21 _ 12_thy_dy - S21_meas; 


F_dx = S11 _22_thy_dx - Sll_meas; 
G_dx = S21 _ 12_thy_dx - S21_meas; 
39 

'/. Compute the partial dervatives 
dFy = (F_dy - F)./delta; 
dGy = (G_dy - G)./delta; 

44 dFx = (F_dx - F)./delta; 

dGx = (G_dx - G)./delta; 


'/. Calculate next permittivity and permeability guess 
xn = x - (1./(dFx.*dGy - dFy.*dGx).*(F.*dGy - dFy.*G)); 

49 yn = y - (1 ./(dFx.*dGy - dFy.*dGx ) .*(-F.*dGx + dFx.*G)); 

% Find the lowest relative approximate error 
epsl = max ( abs (( sqrt (xn) - x))); 

eps2 = max ( abs (( sqrt (yn) - y))); 

54 

eps = max ([epsl eps2]); 

°/ 0 Set the new permittivity and permeability guesses 
x = xn; 

59 y = yn; 
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dy = y + delta; 
dx = x + delta; 

k = k + 1; 

end 

% Pass out the relative permittivites and permeabilities 
epsilon = y; 

69 mu = x; 

Listing A.7: This file uses a variation of the Newton-Raphson 2D root search algo¬ 
rithm. This hie allows for position independence. 

(appendixl / poslndependent.m) 

1 7, Position Independent Analysis 
7, Lt William Gunn 

7 . 

7. The pos Independent () function numerically seeks the complex 
7. permittivity and permeability for the given S-Parameters ... 
regardless of 
6 7. sample position. 


function [epsilon, mu] = posIndependent(Freq,length_Sample ,... 
Sll_meas, S12_meas, S21_meas, S22_meas, x, y,K_csq) 

11 

7. Initialize constants 
tol = 10~(-7) ; 
conv = 1; 

16 k = 0; 

delta = 10'(-7) ; 

dx = x + delta; 

21 dy = y + delta; 


while ((conv > tol) && (k < 100)) 

1 Get S-parameters for standard and derivative functions 
26 [Sll_22_thy, S21_12_thy] = computeSParameter(Freq,. . . 

length_Sample,x,y,K_csq); 

[S1l_22_thy_dy,S21_12_thy_dy] = computeSParameter(Freq,.. . 
length_Sample,x,dy,K_csq); 

[S1l_22_thy_dx,S21_12_thy_dx] = computeSParameter(Freq,.. . 
length_Sample,dx,y,K_csq); 

31 7. Define the standard functions 

F = S21_12_thy.“2 - S21_meas.*S12_meas; 

G = SIl_22_thy.“2 - Sll.meas.*S22_meas; 
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36 


"/. Define the derivative functions 
F_dy = S2 1 _ 12_thy_dy . "2 - S21_meas.*S12_meas; 

G_dy = S11 _22_thy_dy . "2 - SIl_meas.*S22_meas; 

F_dx = S21_12_thy_dx."2 - S21_meas.*S12_meas; 

G_dx = S11 _22_thy_dx . ~2 - S11_meas.*S22_meas; 

% Compute the partial dervatives 
dFy = (F_dy - F)./delta; 
dGy = (G_dy - G)./delta; 

46 dFx = (F_dx - F)./delta; 

dGx = (G_dx - G)./delta; 

"/, Calculate next permittivity and permeability guess 
xn = x - (1./(dFx.*dGy - dFy.*dGx).*(F.*dGy - dFy.*G)); 

51 yn = y - (1 ./(dFx.*dGy - dFy.*dGx ) .*(-F.*dGx + dFx.*G)); 

1 Find the maximum difference between the new and original 
roots 

convl = max(abs((xn - x))); 
conv2 = max(abs((yn - y))); 

56 

conv = max ( [convl conv2]); 

'/. Set the new permittivity and permeability guesses 
x = xn; 

61 y = yn; 

1 Calculate the new deltas 
dy = y + delta; 
dx = x + delta; 

66 

k = k + 1; 

end 

°/ 0 Pass out the relative permittivites and permeabilities 
71 epsilon = y; 
mu = x; 
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